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Due the chiral nature of the Dirac equation, that governs the electron dynamics in graphene, the 
overlying of an electrical superlattice (SL) can open new Dirac points on the Fermi-surface of the 
energy spectrum. These lead to novel low-excitation physical phenomena. A typical such system is 
neutral graphene with a symmetrical unidirectional SL. We show here that in smooth SLs, a good 
mathematical description for particles can be achieved by the semiclassical approximation. Due to 
the one-dimensional nature of the unidirectional potential, a wavefunction description leads to a 
generalized Bohr-Sommerfeld quantization condition for the energy eigenvalues. In order to pave 
the way for the application of semiclassical methods to general two dimensional SLs, we compare 
these energy eigenvalues with those obtained from numerical calculations, and with the results from 
a semiclassical Gutzwiller trace formula via the beam-splitting technique. Finally, we calculate 
conductivities in general point-symmetric unidirectional SLs with one electron and one hole region 
in the fundamental cell showing only Klein scattering of the semiclassical wavefunctions. 

PACS numbers: 03.65.Sq, 72.80.Vp, 73.21. Cd, 73.22.Pr 



I. INTRODUCTION 

Suspended graphene samples exhibit high electron mo- 
bilities, where ballistic transport is seen for samples up 
to the micron length PHH . Within the tight-binding ap- 
proximation, the graphene system shows two inequiva- 
lent momentum energy knots in the Brioullin zone at 
low energies. They are located at momenta K and K'. 
An effective low-energy description around these points 
is given by a massless Dirac equation. Electrons close to 
these points are related to each other by time-inversion 
symmetry The effective quasi- relativistic Hamilto- 
nian is then given by 



H = Hvp —<J T d- 



-Uydy 



V(r) 



(1) 



for electrons near the K point. Here vf is the electron 
velocity in graphene and V(r) denotes an external poten- 
tial. 

In the energy spectrum of the associated Schrodinger 
equation, minigaps arc opened by the application of an 
overlying superlattice potential. New bands arise and 
two of them may touch each other at certain momenta, 
showing up in new Dirac points. These points are classi- 
fied by their local similarity of the energy spectrum with 
the spectrum of the massless Dirac equation. Beside this, 
they show locally a chiral behavior in the pseudospin ex- 
pectation value ((a x ),(a y )) as a function of the Bloch 
momentum Q. 

This was first claimed theoretically by Park et al. @,0] , 
and new Dirac points were in fact found experimentally 
quite recently for graphene on a hexagonal boron nitride 
substrate [§[. One finds new Dirac points at momenta 
G m /2, where G m is a reciprocal lattice wavevector of the 
SL. Their energies are En = ±hvp\G m \/2 [f|. Due to 
their nonzero energies, these new Dirac points cannot be 



observed experimentally in low-energy excitation experi- 
ments on neutral graphene. 

Later on, additional new Dirac points were found in 
theoretical analyses. They were all located at zero en- 
ergy. The calculations were done for graphene with a 
superimposed unidirection electrical superlattice poten- 
tial 0,01 . Actually, these had already appeared earlier in 
the literature within the framework of an unidirectional 
SL on a nanotube [To| . For the most simple representa- 
tion of a unidirectional SL step-potential V(x) = V\{x), 
where x{ x ) = sg[sin(27rx/d)], and sg[x] denotes the sign 
of x, the lowest energy band is shown in Fig. 1. The 
curves were obtained from a precise numerical diagonal- 
ization. The full lowest band energy spectrum shows a 
mirror symmetry at the p Xl Py axes. 

The energy spectrum close to the Dirac points is given 
byPEl 



Wy 



(2) 



p 2 y ) 1 / 2 d/2h, r = sin[a ]e Mo /ao, 



with a = {[V/vf] 2 
&o = ao/V where V = Vd/2hvp. The Bloch momenta 
in x-dircction lie in the Brillouin zone —n/d < p x /H < 
iv/d. The parameter s distinguishes the conduction band 
(s = 1) from the valence band (s = —1). From this 
we deduce that the Fermi- velocities v = Hde s /dp at the 
central valley Dirac point is in general anisotropic [l3j |. 
The new Dirac points are located at momenta p x = and 
p y d/2h = ±(V 2 — (to) 2 ) 1 ' 2 with integer n <E N where 
p y has real values. Furthermore we obtain for momenta 
beyond the new Dirac points \p y \ 3> V/vp the energy 
e s ~ sv F \p y \ [l2|. 

These new zero-energy Dirac points differ from the 
above mentioned points of Park et al. (fjj. They are 
not located at momenta at certain fractions of the recip- 
rocal lattice. It is well known that these momenta are 
part of the region where the SL minibands are formed. 



2 





1 /I 

1 .4 




1.2 


PL 


1.0 












0.8 


MJ 






0.6 




0.4 




0.2 




0.0 



v=o 



7T / 
I 

I 

I 

I 



2tt; 



3tt 



i 
; 

K K 

I J 

If 



4tt 



\// \ 

■n ; v... 



/ 



\ 



10 



12 



14 



p y d/2h 



FIG. 1. Lowest energy band ed/2hvp for the Bloch momen- 
tum p x — and certain SL potentials = 1^x( a; ) as func- 
tion of the dimensionless transversal momentum p y d/2h. 



The Dirac points are then touching points of two mini- 
bands. In the case of the zero-energy Dirac points for 
unidirectional SLs, the pristine Dirac cones are deformed 
strongly due to electrical potential such that the electron 
and valence bands are touched. 

These new Dirac points are located at zero energy, and 
for that reason they posses a number of new interesting 
transport properties 

H El EMI- By application of 
an additional magnetic field, one finds for example new 
Quantum- Hall plateaus In disordered SLs, there may 
also exist interesting localization phenomena [19j . 

A general understanding of the energy spectrum, es- 
pecially the location of new Dirac points, for general two 
dimensional non-unidirectional potentials is still missing. 
Most interesting is thereby the low-energy sector of the 
energy spectrum in neutral graphene. By taking into ac- 
count the fact that the zero-energy Dirac points show 
up only at large SL potentials one may use semiclassical 
methods to determine the lowest energy bands for general 
smooth SLs. For a justification we should demand that 
for unidirectional potentials the semiclassical condition 



Hvf 



\(E-V)(V')\ 



y/{E-v)*-(v FPy y< 



< 1 



(3) 



should be fulfilled except of a few penetration points 
where 



y/[E-V(x p )]*-(v FPy y< 



0. 



(4) 



Within classical mechanics , these points correspond to 
turning points. Due to the chiral nature of ([1]), this is 
no longer true. For example at transverse momentum 
p y = 0, the transmission probability through the pene- 
tration points is unity. Here the electron transforms from 
a particle state (electron-like) to a hole state (positron- 
like )or vice versa. This is the analog of Klein's paradox 
[20l422| in relativistic quantum mechanics. In this context 



it was shown later by Sauter [23|, that the transmission 
probability is decreasing for non-zero momenta p y where 
it approaches zero for VFPy/ftV'(x p ) 7S> 1. 

In order to prepare the semiclassical approach for 
smooth SLs we shall first discuss in Sect. II the semiclassi- 
cal wavefunctions for the problem. Then we derive trans- 
mission and reflection coefficients for electrons or holes 
carrying out Klein's scattering analysis through classical 
forbidden regions between two penetration points. We 
apply our results to the simplest unidirection SL with 
one electron and one hole region in the fundamental cell 
showing only Klein scattering. Our scmiclassic results 
for the lowest energy bands compare well with those ob- 
tained numerically. 

Furthermore, we address in Sect. II the question 
whether one can construct within the semiclassical ap- 
proximation SLs which show the ability to focuse elec- 
tron beams. In order to see in which way a semiclassi- 
cal approach could work also beyond the unidirectional 
SL case we consider in Sect. Ill the generalization of 
the Gutzwiller trace formula that includes also the beam 
splitting phenomenon in the calculation of the semiclassi- 
cal density of states. When calculated from small- length 
orbits, this density of states will permit us to reconstruct 
the lowest lying energy band. Finally we shall calculate 
in Sect. IV semiclassical conductivity formulas for the 
smooth unidirectional SLs and compare our results with 
other known calculations in the literature. 



II. THE SEMICLASSICAL LOWEST-BAND 
ENERGY SPECTRUM 



In the following we formulate the semiclassical ap- 
proach to the quasi-rclativistic Dirac equation of elec- 
trons in a unidirectional SL JI}. The energy spectrum 
will show mirror symmetry with respect to the transver- 
sal momentum p y at the p y = axis. In the first subsec- 
tion, the semiclassical solution of the eigenvalue problem 
will be obtained via the Bohr-Sommerfeld quantization 
condition for non-relativistic electrons. We get very good 
agreement for the lowest band energy spectrum with nu- 
merical results for the SL-dcformed sinus potential. In 
the second subsection, we consider an interesting counter 
example where the semiclassical approach fails. 



Generalized Bohr-Sommerfeld formalism 



Starting point is the semiclassical wavefunction so- 
lution for the Hamiltonian (QJ. This was formerly 
done in the case of the relativistic Dirac-equation in 
Refs. [U,[2^]. One uses the semiclassical Ansatz ^f s (x) = 
~Yj n K nii & n e lStyX ^ h in which ^f n , S(x) are independent of 



3 



h. Up to the order h° we obtain 



\E-V{x)\(s 



\/pI+p 
1 



v , e iS(x)/h+i!p(x) e ipyy/h 



with 



Px{x) 
S(x) = 
0(x) = 



U F 

dx'p x (x') 



-?v r dx > 1 9 X ,(E-V(x')) 



Px {x<) (E-V(x')) 



(5) 

(6) 
(7) 
(8) 



Here s — sg[E — V(x)}, and S(x) is the classical eikonal 
action of the particle. For neutral graphene, the states 
with s = 1 are particle, and those with s = — 1 are hole- 
like. The symbols p x , p y denote the momenta of the par- 
ticle or hole in x, y-direction. Note that for = and 

neglecting the vector part (s(p x — ip y )/ ^P x + P y > 1) T m 

([5]), the wave function *& s (x) is the semiclassical solu- 
tion of the massless quasi-relativistic Klein-Gordon wave 
equation. This means that 4>(x) is a phase correction 
factor due to the chiral nature of the quasi-relativistic 
Dirac equation ([1]). This phase factor has, of course, di- 
rect consequences on the semiclassical Bohr-Sommerfeld 
quantization condition (25l [26j . Without proof we note 
that by taking into account also a homogeneous magnetic 
field this factor cancels exactly the Maslov index |27[ of 
the turning points such that the Landau level energy lad- 
der starts at zero energy [28 1. 

We consider in this paper the simplest case of small en- 
ergies \E\ -C max[|F(x)|] with the scattering in a smooth 
SL is mainly based on the so-called Klein tunneling for 
Py < max[|^(x)|]/i'p'. We shall discuss the situation in 
general also for larger p y at the end of Sect. IIA. 

We assume in the following that the scattering process 
has incident particles from the left side of the tunnel re- 
gion with a positive velocity where E — V{x) > 0. The 
particle tunnels from the left penetration point x p l into a 
classical forbidden region between the penetration points 
x p l and x p r, where \J[E — V{x)] 2 — (vFP y ) 2 is imagi- 
nary. Beyond the penetration point x p r on the right 
hand side of the tunnel region, i t will reach the classi- 
cal allowed hole region, where ^J[E — V(x)] 2 — (vFPy) 2 
is again real but now for E — V(x) < 0. Note that be- 
side Klein tunneling, there are also conventional tunnel- 
ing processes, for example where a particle (hole) tun- 
nels through the full hole (particle) region at imaginary 
y/[E — U(ir)] 2 — (vppy) 2 , i.e., where the particle (hole) 
does not change its signature s. 

By comparing ([5])-([8]) with the definition of the pene- 
tration points @, we obtain a singular behavior of the 
semiclassical wavefunction at these points. This means 
that there is still the freedom to linearly combine the ba- 
sis of semiclassical wavefunction solutions in ([5]) consist- 
ing of left- and right-moving particles or holes with some 



not yet determined possible numerical prefactors in ev- 
ery nonsingular potential sector. This freedom has to be 
fixed by further physical arguments. As in the quasi non- 
relativistic case, we achieve this by matching the wave- 
function ([5]) to the x — > ±oo asymptotics of the exact 
solution of ([1]) for the linear potential V(x) = P(x — x p ) 
where P ss dx p V(x p ) = V'(x p ), with x p = (x pR + x pL )/2. 
In the following we assume that the SL is smooth in 
the classically forbidden region meaning that V'(x) = P 
changes little between the left and right penetration 
points x p L and x p r where P > 0. In order to solve 
this linear potential scattering problem we first use the 
Ansatz *&(x) = hvp[jcr x d x + <r y p y /h — Px](l, l) T ip(x). 
This leads in fact to a Klein-Gordon-like differential equa- 
tion for the wavefunction ip{x) in ([1]). It is given in 
rescalcd coordinates by 



1 ~9 

2 +p y 



(p(x') = 



(9) 



where the dimensionless transversal momentum square 



py is given by 



Pi 



VFPy 

2h\P\ 



(10) 



and x' = (2P/Hvf) 1 ^ 2 x. It can be solved with the help 
of special functions (29|. By comparing the asymptotics 
of this solution for x — > ±oo with the semiclassical wave- 
function ([5]), wc obtain the transmission and reflection 
coefficients (fTTj) . (fT2"j) of the scattering of an electron in- 
cident from the left into the potential V(x) = P(x — x p ). 
After a lengthy but straight forward calculation, the fol- 
lowing reflection i? c h and transmission T^ c h coefficients 



arc found 



i? ch = e^Vl-e- 2 ^ 



with 



0(F?) = -t/4 + arg[L(^)] + p 2 y - f y \n(p 2 y ) 



(11) 
(12) 

(13) 

Note that a similar calculation for smooth gra phene np 
and npn junctions was carried out in Refs. [30l [3l| . We 
show in Fig. 2 the functions \T\ = |7* oh |(p2), \R\ = 

|^eh|(Py) and 3tt/A + ^(p 2 ). Wc obtain that ^ is an in- 
creasing function of p 2 with limiting values i9(0) = — 3ir/4 
and linip^oo d(p 2 y ) = -tt/2. 

The transmission and reflection coefficient for a right- 
incident hole but still with V'(x p ) > can then deter- 
mined from (TTT|) and (fT2|) . using the invariance of (fT]) 
under the transformation 'S'(.t) — > a z ^*(x) for fixed p y . 

This leads to V oh = 7^ oh and ^ oh = ^* h . The trans- 
mission and reflection coefficients for a potential with 
V(x p ) < can be read off from the above coefficients 
by using the substitution p y — > — p y in the correspond- 
ing coefficients. This leads to T~\ lc (p y ) = 7^ c h( — Py) and 
^ho(Py) = R e h(-p y )- 
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Py 

FIG. 2. Transmission function \T\, reflection function \R\. 
and 3^/4 + 1? as a function of fiy . 



oo, finally, we obtain from (|lip . 



For H — > or p y 

(p| and (I3J) that 7* ch -> and ^ e h = e~ i7r / 2 . Thus 
in this limit, we obtain the same reflection and transmis- 
sion coefficient as for the reflection of a non-relativistic 
particle at a smooth potential barrier (32j . 

The matching procedure used above for the determi- 
nation of the reflection and transmission coefficients re- 
quires that the potential changes little between the pen- 
etration points in the classical forbidden region. This 
has to be fulfilled even in the vicinity of the penetration 
points. We find from Eq. that V(x) should be almost 
constant where 



\x - xJ 5, 



] 

p 



1 



-1/3 

Py J 



(14) 



It is clear that a step- like SL does not fulfil this condition. 
Note that in the following subsequent numerical calcula- 
tions, we shall use P~\V(x p r) — V(x p l)\/ (x p r — x p l) for 
fixing p\. 

Next we calculate the energy spectrum for an unidirec- 
tional supcrlatticc potential with two penetration points 
in the fundamental cell. This configuration is shown in 
Fig. 3. In order to calculate the eigenvalue spectrum, we 
have used a transfer matrix method. Note for example 
that the transfer matrix across the i-th penetration point 
is 



M, = 



With the definitions 




(15) 



Si(x) = (-1) 



and 



KM 



3 <tK*) o \ 
e^T^) / ' 



(17) 




(^2,^2 



FIG. 3. We show the deformed sinus potential ()19p for V — 4 
and a = 0, 0.1, 0.5, 100, respectively. The various penetration 
points are shown for the non-deformed sinus potential a — 
in the case that E = and p v vf/V = 1. 



the energy spectrum is given by det[A — e lp * d l h E\ = 
where A = M 2 N 2 (x v ->.>>)M 1 N 1 (x r i 9 ) ■ The various inter- 
section points x pij are illustrated for a sinus potential 
and E = in Fig. 3. We obtain from ©, that the phase 
factors 4>(x) in the transfer matrix N_ { are cancelled. By 
using once more the arguments below (|13j) we obtain 

/S , i + 5 2 \ D n D {S1—S2 
cos I ; 1 — \Rx cos [ : har. 



g^-arg^] 



ITillTalcos 



Pxd 



(18) 



where Si = Si(x p i2). In order to obtain a particle-hole 
symmetry in the spectrum which is the requirement to 
find new Dirac points for E = 0, we demand the point 
symmetry of the SL-potcntial meaning Si = —S2 for 
E = 0. By using the fact that (1 - | J2 X 1 1 J2 2 | )/ |Ti 1 1^2 1 > 1 
for 7^ I-R2I, we obtain that Eq. (fT8|) can not be ful- 
filled for E = in the case of a potential which has no 
additional mirror symmetry with respect to the transver- 
sal momentum p y . This means that semiclassically we do 
not find any additional Dirac-point except that for pris- 
tine graphene at p x = p y = for asymmetric potentials 

where We can see this numerically by cal- 

culating for the deformed sinus-potential 



V(x) = Vsin 



2tt(x - d/2) 1 + a(x - d/2) 2 ) 



l + a{d/2f 



(19) 



defined for < x < d, the lowest energy band by a 
numerical diagonalization, and compare the results with 
the semiclassical ones (fT8|) . In Fig. 3 we show various 
deformed sinus potentials for which we did that. 

We show in the left panel of Fig. 4 the lowest energy 
band for p x = by using the exact numerical diago- 
nalization method for the various sinus potentials fj 19|) . 
In the right panel the corresponding semiclassical result 
(|18p is shown. The same is shown in Fig. 5 for various de- 
formed sinus-potentials. The plots are characterized by 
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sults (Tig)) . 



2 4 6 



10 120 2 4 6 8 10 12 



2.0 

ft, 

9 

* 15 



0.5 



0.0 



ex 


> { 

. / 

; / 


SC 






10 ?''' / \\, 


^, iii 

\ • i ■ 




1 

; i 


1 '. \ A 


/ /o.5 r 


\ /// 

i . ; i / 




i 


i \ \ • ,' 
i\ \ \ • i 


//"X i • 

\ / 


> \ . / . 

v ' ,' 1 i 

, 1 V ' ; 

', \ ' ' 






i \ \ i / / 
/ • \ \ ' / 
/ \ \ V / ■ 


/' a=fl.l^' 

/ 


'. >,'',/ / 
'. V » ' 

■v\/ 


t 




/ ' \ 'V ' 

/ '-\ t\ ! 
w »/ 



2.0 



0.5 



0.0 



2 4 6 8 10 120 2 4 6 8 10 12 

p s d/2/S p,rf/2S 

FIG. 5. Lowest energy band at p x — for the deformed sinus- 
potential (fT9|) with V" = 47T and a = 0.1,0.5, 100. Left panel 
shows the numerical, right panel the semiclassical results (|18l) . 



V(x) = V(x)d/2hvF since the energy spectrum (up to 
a simple rescaling of momentum and energy) as well as 
the conductivities depend mainly on this dimensionless 
potential. From both figures we obtain an almost perfect 
agreement between numerical diagonalization results and 
semiclassical lowest energy band. 

As was already discussed below ((SJ), the Klein- 
scattering process is dominant for small energies E and 
momenta p y over other scattering processes for smooth 
SLs. Let us elaborate this point further. First we can 
show using semiclassical methods similar to those ap- 
plied to the step-like case in Sect. I, that |e s | > vf\p v \ 
at large momenta where |p y w_F| m&x[\V(x)\]. On the 
other hand for small momenta where |_Pj/^f| is smaller 
than the absolute value of possible local minima (max- 
ima) of V(x) in the case of particles (holes), we obtain 
for small energies \E\ <C Ip^^fI that mainly Klein scat- 
tering processes are active. For momenta p y between 
these two extrema, also conventional scattering are rele- 
vant. To avoid them at low energies we must take into 
account (|14|) . and demand that the SL potential V{x) has 
no local minima for particles and maxima for holes and 
that the local minima and maxima are of similar abso- 
lute potential value. Furthermore we must demand that 
V'(x) « const within the local minima and maxima. The 
smooth forms of the symmetric two-step potential are in 



the class of potentials fulfilling these requirements. 

We point out that these requirements are not necessary 
but sufficient to determine the whole low energy region 
of the lowest energy band for a given SL potential within 
the semiclassical method discussed in this paper. The 
reason that these requirements are not necessary lies in 
the fact that the type of scattering depends strongly on 
the energy that the particle or hole has. The above re- 
quirements held under the assumption that \E\ <?C \p v vf\, 
and this must not be fulfilled for certain momentum val- 
ues p y . 

In this context let us mention that one can set up an 
SL potential V(x) which shows a new Dirac point at zero 
energy, where conventional scattering takes place. Here 
we use that for potentials where the tunnelling regions 
and the effective tunneling potential is large enough, the 
particles and holes are fully backscattered at the penetra- 
tion points. Then we can consider the classical allowed 
regions as independent. In the case that the semiclassical 
condition S/h — n/2 = n(n — 1) with n € N is fulfilled in 
every classical allowed region, where S is the eikonal in 
the region, we find a new Dirac point in the spectrum. 

In the field of quantum chaos of non-relativistic sys- 
tems, the energy spectrum is very important to charac- 
terize it. This is seen for example in the nearest neigh- 
bour level statistics via a non-Poissonic behavior [331 ] . 
Such a behavior can be also seen in the energy spectra 
of random matrices. The main difference between the 
quasi-relativistic SL systems we are looking at and the 
quantum chaos in non-relativistic systems lies in the ener- 
gise of a few lowest energy bands labelled by the Brioullin 
zone wavevectors (p x ,p y ) with \E\ < max[\V(x)\]. These 
are the relevant bands governing the most interesting low 
excitations in graphene at half filling especially, where we 
find new Dirac points. At larger energies, no backscat- 
tcring is seen semiclassically in the SL system. This has 
the consequence that minibands together with the cor- 
responding minigaps are not shown up at high absolute 
energies in semiclassical approximations. 



B. Constructing SLs for focusing electron beams 

In the following, we restrict ourselves now to unidi- 
rcction SLs with a mirror symmetry and an additional 
point symmetry at the origin similar to the sinus poten- 
tial discussed above. As was shown in the last section, 
this requirement is necessary to find new Dirac points on 
the E = axis. Then we obtain from (fl~8f for E = 0, 

S = Si = —S2 and arg[i?] = arg[i^i] = — arg[i^ 2 ], that 
the momentum p y = p™ of the new Dirac points is deter- 
mined by 



S 



arg 



(20) 
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1.5 

p y d/2H 

FIG. 6. The (black) solid curve (sc) is the semiclassically 
calculated energy spectrum for the potential Vf at p x = 
which is shown as the black solid curve in the inset. The other 
curves (ex) in the main panel are the lowest energy spectra 
calculated by an exact numerical diagonalization calculation 
for the potential Vf and its variations shown in the inset. 



where n£fj. The number of new Dirac points is then 
given by 



1 



IS 

7T H 



arg 



1\ 



(21) 



where we must set p y = in (|2Tj) . We have used the 
abbreviation [x] as the largest integer number smaller 
than x. 

As mentioned above and can be deduced from for 
the unidirectional step- like SL, electrons with a momen- 
tum near the central Dirac point are focused strongly in 
the direction of the SL wavevector, i.e., v x ^> v y , espe- 
cially where a new Dirac point emerges. It was mentioned 
in Refs. [IH that this phenomenon could have technical 
applications for a strong focusing of electron beams in 
graphene. Of course, a true focusing of an electron beam 
demand the additional requirement that v y = 0, in the 
vicinity of some momentum value and not only exact on 
that value. Such energy dispersions were in fact found in 
photonic crystals [H, |35| . By solving (|20|) in a nontrivial 
momentum region we are now able to construct within 
the semiclassical approximation potentials showing also 
such a behavior. For doing this we restrict ourselve in 
the high potential region to the following sort of poten- 
tials V sin[27r(.T - d/2 - 5d)/d] for d/2 < x x < x < 3d/4. 
Here we can restrict ourselve in the discussion to the 
positive branch of the potential, i.e. d/2 < x < 3<i/4, 
being sufficient due its symmetry. The value x\ is given 
by the condition that (|2T)|) is fulfilled for the momentum 
p y = V(xi)/vp where we restrict ourselve in the follow- 
ing to n = 1. 

We now determine the potential Vf(x) for d/2 < x < 
xi by solving (|20|) itcratively. Here we determine Sd and 
V so that Vf(d/2) « and further that momentum value 
of the first Dirac point p y is maximal. With these require- 
ment we obtain 2irSd/d = 0.437 and V = 4.58. 



This then leads within the semiclassical approximation 
to the fact that the first side-valley Dirac point p y and 
the central Dirac point is connected by a flat energy dis- 
persion curve with zero energy. We show in the inset 
in Fig. 6 by the black solid curve the potential Vf(x) 
obtained in this way. The (red) dotted and the (blue) 
dashed potential curves are variations of Vt being dis- 
tinct at x < x\ -values. The black solid curve denoted 
with sc in the main panel in Fig. 6 shows then the semi- 
classically calculated lowest energy spectrum by using 
(|18[). We obtain in fact a flat energy curve around the 
central Dirac point. The other energy curves shown in 
the Figure are calculated by using the exact numerical di- 
agonalization method. The various exact diagonalization 
curves in Fig. 6 correspond then to the potential varia- 
tions shown in the inset. We obtain from the (black) 
solid energy curve that the flatness seen in the semiclas- 
sical approximation is in fact vanished within the nu- 
merical diagonalization calculation. Note that this even 
holds when going to a high basis number in the exact nu- 
merical diagonalization calculation. This shows that the 
semiclassical approximation fails here for the constructed 
potential Vf at least for the effect in showing the flatness 
of the energy spectrum close to the central Dirac point . 
The reason for this failure comes presumably from the 
fact that at the penetration point x P 2i = x\ the condi- 
tion (fT4"|) is no longer fulfilled. We even obtain from Fig. 6 
that the energy curves of the potential variations of Vf 
shown in the inset, do not vary much around the central 
Dirac point. We consider this as a hint that presumably 
the whole program in finding an SL potential having a 
flat piece in the energy spectrum, with one electron and 
one hole region in the fundamental cell, seems doomed to 
failure. Here we take also into account that Vf even gets 
more shallow for larger x\p y in the range d/2 < x < X\. 
Note also that we carried out further numerical calcu- 
lations with variations of the SL potential which were 
also not successfull. It was shown in Ref. [T3 that such a 
program can be successful when considering more com- 
plicated SLs. In that paper it is shown that an SL with 
one electron and one hole region and an additional small 
modulation of the potential strengths over many funda- 
mental cells of the SL can lead to energy spectra with a 
flat behavior around the Dirac points. 



III. SEMICLASSICAL DENSITY OF STATES 

The generalization of the above results to general two- 
dimensional SLs via a semiclassical wavefunction solution 
of Eq. ([1]) is not possible. In the case of non-relativistic 
quantum mechanical systems this can be carried out only 
for integrable systems j36f . This result is modified in 
relativistic systems mainly due to the existence of the 
additional phase factor 4>{x) in ([5]) [25|, [26J. One way out 
of this dilemma is by calculating the density of states 
semiclassically with a formalism developed by Gutzwillcr 
[3?J ■ The eigenvalue spectrum is then determined from 
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the calculated density of states. 

For the unidirectional SL with one electron and one 
hole regionper fundamental cell we obtain for the density 
of states Hi, [H] p(E) = p{E) + Ap(E). Here p is the 
average density of states. It is given by 



1 f d 
p(E) m — / dxRe 



\E-V(x)\ 



2 /(E-V(x))i 



P7, 



(22) 



The fluctuating part is given within scmiclassical approx- 
imation by 



=i 



Ap{E ) = ±-Re 
nn 

p V—Y 

(23) 

The sum p in (|23p runs over the primitive periodic orbits 
of particles E — V(x) > or holes E — V(x) < 0, respec- 
tively Particles and holes are transformed to each other 
at the penetration points. The configuration space of the 
orbits is given by the fundamental cell of the SL with peri- 
odic (circle) boundary conditions. T p is the required time 
for the particle or hole for passing the primitive orbit. 
T T W stands for = ^iM^Cri^tof^Cri 

Here ~^i{p) (Vi(f))) is the number of transmissions from 
left (right) to right (left) through the potential barrier i 
in the primitive orbit. is the corresponding total 

reflection coefficient. S p is the eikonal of the primitive 
orbit, i.e. S p = n-iSi + n 2 /S 2 where n\ and n 2 are the 
number of transitions of the particle regions E—V(x) > 
and hole regions E — V(x) < 0. j(p) is the winding num- 
ber of the primitive orbit on the circle representing the 
fundamental cell. 

In order to derive Eq. ([23jl. we used the ray split- 
ting generalization of Gutzwiller's trace formula first dis- 
cussed in Ref . H^. There it was shown that a ray splitting 
boundary in an integrable system can cause additional 
sign of chaos in the energy spectrum. One of the simplest 
systems with ray splitting, is that of a non-relativistic 
electron in an infinite one-dimensional square well with 
a discontinuous step inside the well [40l - f42l ]. This sys- 
tem but also our system represented by the density of 
states (|23|) can be discussed in the background of quan- 
tum graphs [43| . With the help of the methods used in 
this reference one can directly show the connection be- 
tween the density of states (|23|) and the corresponding 
energy spectrum represented by equation (|18|l . 

By using that the absolute value of the particle veloc- 



ity is given by Vf^J(E — V(x)) 2 — p 2 v 2 F /\E — V{x)\ one 

can easily determine T p for every primitive orbit by in- 
tegrating the inverse velocity over the orbit. It is well 
known [36[ that a renumbering of the summands in (|23|) 
can lead to divergent subseries. A well behaved approx- 
imation should be reached by sorting the terms in (|23|) 
with respect to their maximal orbit length l max . This 



equal to 2md. We show in the upper row in Fig. 7 Ap 
f° r 'max = 2d , 4d and 8d for the sinus potential (fT9|) 
with V = 27r, a = and p x = 0. The various panels are 
calculated for p y d/h values where the lowest energy band 
(cf. Fig. 4) has its two Dirac points (left and right panel) 
and further where the band has its local maximum (mid 
panel). By comparing these panels with the correspond- 
ing energy spectrum in Fig. 4 we obtain in fact that at 
least approximately the energy values of the lowest en- 
ergy band shows up by the smallest energy maximum 
in Ap. This happens when we even take into account 
only small l max orbit lengths. We note that the higher 
energy maxima in Ap in Fig. 7 correspond to higher en- 
ergy bands. Next we try to reproduce from the lowest 
energy maximum in Ap, the lowest energy band for var- 
ious (deformed) sinus potentials and Z max = 2d, 4d and 
8c?. We compare our result in Fig. 7 with the semiclassi- 
cally calculated lowest energy band by using (fT8|) ((black) 
straight curves) . We obtain from the figure that the new 
Dirac points even show up for small Z max by a plateau 
at zero energy where its extension is rapidly decreasing 
for higher Z max - values. Note that the maximum criterium 
used here for determing the spectrum from Ap is distinct 
from the common approaches in the background of de- 
terming the full energy spectrum for systems in the field 
of quantum chaos. There one commonly uses the fact 
that by integrating the full density of states between two 
non-degenerate energy levels gives the value one. Since 
we are only interested in the lowest energy level and fur- 
ther since the lowest energy band and the first excited 
energy band is well separated such an approach is not 
necessary here. 



IV. CONDUCTIVITIES 

Next, we calculate the conductivities parallel and or- 
thogonal to the SL wavevector by using the semiclas- 
sical wavefunction ([5]) and energy dispersion (|18|) for 
SLs with a point symmetry at zero energy for half- 
filling. The lowest band eigenvalue spectrum is given 
by (Tj"8]) which was effectively calculated from the eigen- 
values of the matrix A. By using (jT5|) we obtain the 
following energy dispersion around the Dirac points, i.e. 
\e a \dI(-l,0)/2v F S 1, 



2hv f 



d J(-1,0) 



|T| 2 sin 2 



2 (S + a r g [l] 



Pxd 
2h 



-(l-|i?| 2 -|T| 2 ) 



(24) 

1/2 



means that for I 



2md with to S N we have to take 



where s = 1 for the conduction band and s = — 1 for the 
valence band. Here we use the abbreviation arg[it] = 

(argf^h - arg[it 2 ])/2 and \R\ = W^lffil, \T\ = 



into account in (|53J) all orbits with lengths less than or J |3*i||7* 2 | and further we will use ^ = |i?|e larg ^] . The 
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FIG. 7. Upper panels show the dimensionless density of states 2AhvFp/d as a function of the dimensionless energy Ed/2hvF for 
various maximal orbit lengths Z max at p x — 0. The SL potential is given by the non-deformed sinus potential (|19[) with V — 2n 
and a = 0. Left panel is calculated for the central valley where p y d/2h = 0, at transversal momentum value p y d/2h = 2.58 
where we found the maximum of the lowest energy band (middle panel) and at the new side valley Dirac-point momentum 
p y d/2h — 4.39 (right panel). Lowers panels show the energy of the lowest maximum value of the density of states Ap for various 
orbit lengths Z max . We compare these values with the semi-classical spectrum calculated by (|18p (black curve). Left panel is 
calculated for the sinus potential (|19p with V = 2n, a = 0; V — 4n, a = (middle panel) and V = 2n, a — 100 (right panel). 



function I(ni,ri2) is defined by 



I(m,n 2 ) = 



dx 



y 1J y 



1 n n 2 
"y 



\V(x)/v F \ ni+ri2 



(25) 



In the following, we will use the eigenfunctions of 
the matrix A which was defined below (fT8|) . These 
arc given in the vicinity of the Dirac points, i.e. 
\e s \dI(-l,0)/2v F S < 1, by 



/sm 



-sin 



i?| 2 +£|^|T| 2, 



e J f |i?ile l f- 



2f+2arg[^ 
[^1 - l-Role -4 ! + a *s[3.2 



(26) 

The lowest-band eigenfunctions for electrons in the SL 
are then given for < x < d by 



u s (x) « \ M.iKi(x)Q(x - x p u)Q(x pl2 - x) 

+ N 2 (x)M 1 N 1 {x p i 2 )e(x - x p2 i)6(a; P 22 - x^BV/N. 

(27) 

Here JV is a normalization constant. Note that we omit- 
ted here once more semiclassical phase factors ([8]) as pre- 
viously in (|T5)) . We show below that they will in fact 
not contribute to the conductivity within the semiclassi- 
cal approximation. Furthermore we idealized in (|27p the 



whole wavefunction in such a way that in the classical 
forbidden region it is given by zero. We will also justify 
below this assumption 

Next we calculate the dc-response in the SL system. 
This is done in the gauge A = — 6EtQ(t). The con- 
ductivity in the i-th direction in the lowest energy level 
approximation valid for t — > oo is then given by }12l |44| 

' f ^Rele-^U-iMu+Oe+W], (28) 

BZ r 



(2tt) 2 Jn 7 .h 2 



with 



dt' 



dt"T{t"), 



(29) 



t'=0 Jt" = -oo 



and the transition matrix element T = e«"" (Ui|o"i|u_i/ 
The value Ae is given by the energy gap Ae = e\ — e_i 
for an electron with momentum p y . The integral in (|28|) 
is carried over the full Brioullin zone. 

In the following we calculate separately the contribu- 
tion of every energy valley to the momentum integral in 



, i-e., 



£*S(2- 

n=0 



(30) 



For large times one can restrict the p y -integrals of 
Eq. (|2"8"j) to the neighbourhood of the valley center in 
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T where p™ is determined by (|20[) for n > and p° = 
for the central valley. The factor two in (|30jl takes into ac- 
count the mirror symmetry of the energy spectrum with 
respect to p y , such that we may consider only p™ > in 
((301). 

For calculating the conductivity er^, we first have 
to determine the matrix element (ui|crj|u_x). We ap- 
ply the semiclassical approximation by assuming that 
h is so small that we may neglect integrals of the 
form J dxe l2 f dx Si< ^ x ^ h in comparison to the integrals 
/ dxe M I dx Si ( x ^>l h . On similar grounds we may also 
neglect the matrix contributions in the classical forbid- 
den regions. From this argument it becomes even clear 
that the semiclassical phases <p(x) ([8]) will not contribute 
to (ux|<Xj|u_i) since the integrand in ((5J) is inverse pro- 
portional to l/p x (x'). 

By using (j2"7| - (f2l?)) we obtain the following conductiv- 
ities 

- °> (31) 



2 h \v x v y \/v% 



with v x and v y are the electron velocities at the Dirac 
point. By using (|24|) we obtain 



\T 



I(-1,0)' 
2vph 



l(-l,0)d 



(32) 
(33) 



x d Py J \R\ 2 sin 2 f| + arg[i£] J +1(1 - \R\* - \T\*) 



where 



#±(i?) = ±1 - \R\ cos(tf) + \T\ sin(#). (37) 



In the following, we specify further the parameters in 
(|31[) . For the side- valleys n > 0, the value p y in the ex- 
pressions (|31 [) - (|36p is given by the side- valley Dirac point 
momentum governed by (j2"0|) . In this case we obtain for 



v y fl32J) 



2v F 



J(-1,0) 



~/(-l,l) + ^arg[^] 



\R\, (38) 



with 
h 

PyVF 



a PB arg[^] = i[Re[*( i p 2 y )]-log(p 2 )] 



^'(x p21 ) y(x pl2 ) 



x p21 ~ x pl2 \ 

d J 
(39) 



where $ is the digamma function. 
For the central valley n = we have 



/ ttHvf 
d 2 V'{x pl ) 



(40) 



Here a; p j is the i-th intersection point of the SL potential 
and x-axis, i.e. x p \ = d/2 and x P 2 = 0. The p y momen- 
tum value in the expressions (|3ip ~ (|39p is then given by 
p y = 0. The electron velocity in y-direction is then for 
n = given by 



The absolute square of the transition matrix elements are 
given by 



2tt 

: '„ M \N d Fl + N nd Fl | \N d F* + N nd F* \ ' 

2-n (N d Fl + N" d F 4 ) 2 

d-d rrrr^ .". „_,,„,, (34) 



(iV^) 2 



v l \N d F 2 + N nd F^ \\N d F 2 + N nd F3 
with 

AT d = |JQ + 2|i?| 2 )/(-L,0) 

+ Re[fl][/(-l, 2) - 7(-l,0)] + 2Im[^]I(0, 1) 



N nd = —^\R\I(-1,0) 



N d = 4/(0,0), 

N v = [^2 { Re t^] W" 1 . !) + M%(0, 0) J 



(35) 



and 



F± = H + (-d)H^(d)± cos 2 (i?) , F| = (tf ) + cos 2 (i?) , 
i^ = 2cos(tf)# ± (#), F 4 = -2cos(t9)(iJ±(i?) + 1), (36) 



27w, 



d/(-l,0) 



,5' 



3tt 

T 



(9 Py \Ri\)(d Pv \R 2 \) 



(dpjR^-d.JRiiy 



1/2 



(41) 



The only non-zero values in (|35p for pj, = are given 



by iV d = 5 and N d = 4. This leads to O y 
O x = 16tt/25. 

For a step-like SL potential V(x) = V^ir) 
conductivities are given by [l2l . Il5l | 



2 



1 



VV 



e 2 7T 1 

X25 2 ' r ' 1 



= and 
the dc- 

(42) 



with a = nn/V, T„ = (V 2 - (irn) 2 )/V 2 , V = Vd/hv F 2. 
The index n denotes the valleys n = 1, . . . , [V/7r], where 
[x] is the largest integer value smaller than x. Here n = 1 
denotes the outermost valley, and n = [V/tt] the first 
valley beside the central one. For the central valley, we 
have r = sin(V)/V and a = 1. 

We show in Fig. 8 the conductivities a xx {o yy ) in the 
left (right) panel as a function of the potential strength 
V for the non-deformed sinus potential (fT9|) with a = 0. 
We obtain from the figure that for a xx the central val- 
ley contribution cr° to the conductivity is most relevant 
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FIG. 8. Left panel shows the valley contribution a^ x to the 
conductivity in parallel direction to the SL wavevector calcu- 
lated within the semiclassical approximation (|31|) as a func- 
tion of the potential strength. Here we used the non-deformed 
sinus potential ()19|) for a — as the SL. Right panel shows 
Oyy for the same potentials. 
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FIG. 9. Conductivity a xx for a deformed sinus potential of 
the form (119|) as a function of the deformation parameter a 
at various potential strengths V . 



where for a yy the outermost valley a yy contributes max- 
imal. This is in accordance to the case of the step-like 
potential V(x) = Vx( x ) (|42|l . We even obtain from the 
figure that in practice one can neglect the non-dominant 
valleys in the calculation of (|30|) . This is in contrast to 
the step-like case where the non-dominant valley contri- 
butions are much larger. The reason lies in the fact that 
for smooth potentials V{x), erg (f3"Tj) contains exponential 
damping terms as a function of p y via their dependence 
on the transmission coefficient |T|. We obtain from (I3T1) 



IT 



3 ~n 



1/\T\, The exponential vanishing bc- 

^2 



havior of the transmission coefficient \T\ for large p y in 
smooth potentials is caused by the exponentially damp- 
ing of the wavefunction in the classical forbidden region. 
In contrast to this, the transmission coefficient \T\ for a 
step-like potential V(x) is only algebraically vanishing as 
a function of p y . 

One can understand the |T|-behavior of (J3J) also in the 
following heuristic way. The finite quantum conductiv- 
ity in pristine graphene is heuristical conceived by taking 
into account Einstein's law for diffusive scattering. There 
the conductivity is proportional to the density of states 
times the diffusion constant. As in every two dimensional 
system for infinite small scattering the effective diffusion 
constant is infinite. By the same time, in contrast to two 
dimensional metals where the density of states is con- 
stant, it vanishes in graphene at the Dirac point leaving 
then the total conductivity as a constant. By the appli- 
cation of an SL in x-direction, the diffusion in y-direction 
is in first approximation as in graphene, but the density 
of states scales with 1/\T\ (J24J) leading to <r yy ~ l/\T\. 
In contrast to this, the scattering in the ^-direction for 
graphene with a superimposed SL is for |T| -C 1 mainly 
diffusive, with a diffusion constant ~ \T\ 2 . The density of 
states still scales with 1/|T|. Since the density of states 
vanishes at the Dirac point we obtain an extra |T| 2 -term 
in a xx leading to a xx ~ 1^1 3 • 

From (|3"Tj) we obtain the fact that even for small SLs 
where only the central Dirac point is present, a yy is zero. 



This is not true for the orthogonal conductivity a yy of 
the step- like SL system (|42|) . In Ref. [H the conductivity 
a X x for the non-deformed finite length sinus SL-potential 
(fT9f as a function of V was calculated by using a trans- 
fer matrix method for heavily doped graphene [45|, |46| 
leads. We obtain a good quantitative accordance of our 
result in the left panel of Fig. 8 with their curves. We 
consider this as a justification of the semiclassical ap- 
proximation method considered in this paper. In this 
context we should note that in the case of the SL-frce 
pristine graphene system the exact calculation by using 
the transfer matrix method 14a, Wfl and the method used 



here [44j which is justified for undoped leads gives similar 
numerical values. 

Finally we show in Fig. 9, a xx for the deformed sinus- 
potentials (fT9|) as a function of the deformation param- 
eter a for various potential strengths V. We argued in 
Sect. II that in this case only the central Dirac point 
is existent leading to the fact that a yy = within the 
semiclassical approximation. We obtain from Fig. 9 local 
maxima in a xx at certain deformation values a. As can 



be seen from <|31 1) with (|4ip . these deformation parame- 
ters are in the vicinity of deformation parameters where 
([20]) is fulfilled for some value n £ N and p y = where 
arg[^] = -3tt/4. 



V. SUMMARY 

We have analysed the behavior of electrons in elec- 
trical superlattice potentials within the semiclassical ap- 
proximation. We found this description to work well for 
smooth superlattice potentials. We started in Sect. IIA 
by representing the semiclassical wave function of the 
quasi-relativistic Dirac equation of electrons in graphene 
superimposed by an SL. We have derived transmission 
and reflection coefficients for Klein tunneling through a 
classical forbidden region, in which a particle state is con- 
verted to a hole state or vice versa. Within a generalized 
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Bohr-Sommcrfcld formalism, we have derived the eigen- 
value equations for the lowest energy band of an SL with 
one electron and one hole region in the fundamental cell, 
showing only Klein scattering. For electrons in an SL 
of a (deformed) sinus shape, we obtain very good ac- 
cordance of the scmiclassical energy spectrum with the 
spectrum obtained by exact numerical diagonalization. 
Then we try to construct in Sect. IIB potentials having 
an energy plateau at zero energy which could be use- 
ful for focussing electron beams in graphene. We were 
successful with this construction within the semiclassical 
approximation. We have shown that the energy plateaus 
found in the semiclassical approximation do not persist 
when going beyond this approximation. The reason lies 
in the fact that an acceptable smoothness behavior of the 
constructed potential which is necessary for the validity 
of the semiclassical approximation was not given. In or- 
der to path the pave to take into account SLs which are 
not unidirectional we calculated in Sect. Ill the semiclas- 
sical density of states within the generalized Gutzwiller 
trace formula by taking into account the beam-splitting 
extension. Even by considering only small length orbits 
we could reconstruct from the density of states maxima, 
the lowest band energy spectrum. This was carried out 
explicitly for the (deformed) sinus potential SLs. Finally 
we have calculated in Sect. IV longitudinal conductiv- 
ities along and transverse to the SL wavevector within 
the semiclassical approximation. Here we restrict our- 
selve again to the simplest point symmetric SLs with 



one electron and one hole region in the fundamental cell 
where only Klein scattering is important. We obtain a 
good quantitative accordance with conductivity curves 
already existent in the literature for sinus potentials as a 
function of the potential strength. In these calculations 
a transfer matrix method was used in order to calculate 
the conductivity parallel to the SL wavevector. Further- 
more we obtain, as was formerly shown also for step-like 
SLs, that the conductivity along the wavevector of the 
SL is mainly governed by electrons in the central valley 
whereas the orthogonal conductivity is governed by the 
conductivity contribution of the outermost valley. The 
contribution of electrons in the central valley is zero in 
the latter case. In contrast to the step-like SLs, the ne- 
glection of the other non-dominant valleys is exponential 
damped in both cases. This is connected to the fact 
that the transmission coefficients for Klein tunneling in 
smooth potentials in contrast to step-like SLs are expo- 
nentially small as a function of the length of the classi- 
cally forbidden region and transversal momentum. 

As argued in the introduction of this paper, due to 
their spectral and conduction properties, electrons in 
graphene with an overlying SL are interesting systems 
promising many applications, that could open new routes 
in building electronic devices. The ability to construct 
for an energy band with desired conduction properties 
SL potentials which show this behavior, can be very use- 
ful. We have shown that this is in principle possible 
within the semiclassical approximation in our example 
in Sect. IIB. 
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